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We present numerical solution of the chlorine dioxide-iodine-malonic acid 
reaction-diffusion system in two dimensions in a boundary-fed system us- 
ing a realistic model. The bifurcation diagram for the transition from non- 
symmetry breaking structures along boundary feed gradients to transverse 
| symmetry-breaking patterns in a single layer is numerically determined. We 

find this transition to be discontinuous. We make connection with earlier 
results and discuss prospects for future work. 
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I. INTRODUCTION 



In 1952, Alan Turing showed |I| that under certain conditions, reaction and diffusion pro- 
cesses alone could lead to the symmetry-breaking instability of a system from a homogeneous 
state to a stationary patterned state. The Turing instability is characterized by an intrinsic 
wavelength, in contrast to hydrodynamic instabilities, such as Rayleigh-Benard convection 
where the wavelength depends on cell height. For this reason, this instability mechanism has 
particular relevance to pattern formation in biological systems [0,0]. Given the difficulties of 
noninvasive experiments on biological systems, experimental studies of Turing pattern for- 
mation have been carried out primarily on chemical systems. Even so, Turing patterns have 
been generated in the laboratory only recently, specifically in the chlorite-iodide-malonic 
acid (CIMA) chemical reaction-diffusion system 

Although Turing patterns have been extensively studied theoretically and numerically in 
the context of abstract models, the possibility of comparison with controlled and reproducible 
experiments provides motivation for quantitative analyses based on realistic models of these 
systems. Lengyel, Rabai and Epstein (henceforth referred to as LRE) have proposed a real- 
istic model of the simpler chlorine dioxide-iodine-malonic acid (CDIMA) reaction-diffusion 
system [0,0. The chemistry of the CDIMA and CIMA systems are related, and the two are 
similar in terms of their stationary pattern forming and dynamical behavior. The potential 
for both experimental and theoretical work makes the CDIMA system well-suited for the 
study of nonequilibrium pattern formation in general. 

In practice, this has not been fully realized, and unlike in fluid systems [|I]], experimental 
and theoretical efforts in chemical systems have not been closely coupled. First, numerical so- 
lution of react ion- diffusion equations using realistic chemical parameters is computationally 
demanding. In addition, the algebraic complexity of the realistic nonlinear reaction terms 
renders these models unsuitable for analysis by standard analytical tools. Furthermore, un- 
like the case considered originally by Turing and subsequently by others, the experimental 
system is not uniform. A continuous supply of reactants is fed into the reactor from two sep- 
arately nonreactive reservoirs to keep the system far from equilibrium, setting up gradients 
in the concentrations of these reservoir species along the width of the reactor. In an earlier 



work ||11|| , we used the realistic LRE model to investigate the formation of one-dimensional 
stationary structures along the boundary feed gradients and their linear instability to trans- 
verse symmetry-breaking (Turing) patterns. Here, we extend these results by numerically 
solving the fully nonlinear reaction-diffusion equations in two dimensions. 

Kadar et al. |T2j have also numerically investigated two-dimensional Turing patterns 
within the LRE model, however in a closed system (i.e., in the absence of gradients) where 
the patterns are transient by nature. Two-dimensional numerical simulations of Turing pat- 
terns in ramped systems have been performed using popular abstract models, such as the 
Schnakenberg model [IH|,I1|, and the Brusselator model [[TSJ^T^] , as well as the generalized 



Swift-Hohenberg equation ||18|| . These models, which are easier to implement, produce pat- 
terns which possess similar features to those observed in experimental systems. However, 
they do not allow quantitative comparison or prediction of new features at specific parameter 
values of a real chemical system. 

In this article, we present the first numerical solution of the LRE model in a boundary-fed 
system in two dimensions, corresponding to sustained patterns in thin-strip open reactors. 
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We determine the branch of the bifurcation diagram corresponding to the transition to 
stripes in this system, a result which can be directly investigated in experiments based on 
the CDIMA system. In Sec. P], we present the model and parameters used. In Sec. PTT| , we 
describe our numerical method. We present our results in Sec. [TV|, and give conclusions and 
prospects for future work in Sec. |V[ 

II. CHEMICAL MODEL 

The LRE model of the CDIMA react ion- diffusion system has been discussed in detail 



elsewhere |T2| , p!T|JT9|] . The resulting governing partial differential equations are given below: 

d[MA] 



—r\ + -DmaV 2 [MA], (1) 

-ri + \r 2 + 2r 3 - r 4 + A 2 V 2 [I 2 ], (2) 

-r 2 + £ C io 2 V 2 [C10 2 ], (3) 

n _ r2 _ 4r 3 - r4 + D T - V 2 [I - ] , (4) 
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The nonlinear reaction terms are given by: 

MMA][I 2 ] 

ri " k lb + [l 2 ) ' (8) 

r 2 = fc 2 [C10 2 ][n], (9) 

ni = Wao a ir]W + M ™ a g ri , do) 

r 4 = fc + [S]p[ 2 ][I-]-A;_[SI 3 -]. (11) 

The rate and diffusion constants used in the numerical calculations here are taken from 
Refs. ||12| , |l9"|p0| and are given in Table |. The left/right reservoir species are malonic acid 
(MA) /chlorine dioxide (C10 2 ) and iodine (I 2 ), respectively. As these species diffuse through 
and react in the gel reactor, the dynamical iodide (I - ) and chlorite (C10 2 ~) species are 
produced. The starch triiodide complex Sl3~ is the experimentally observed species. 

The experimental CIMA reaction is similar to the CDIMA reaction in terms of dynamics 
and stationary pattern formation |i~9|,pl|,p2[ . However, a quantitative model of the CIMA 
reaction does not exist. In particular, the role of chlorine-dioxide in this reaction is not 



well understood [19|. Furthermore, it has been pointed out that in the experimental CIMA 



system, use of chlorite and iodide along with acid as reservoir species could lead to reactive 



reservoirs [13 . In the CDIMA system, chlorine dioxide is non-reactive with both iodine and 
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malonic acid, allowing for well controlled boundary concentrations of these species. It is 
known from batch experiments and simulations of the CDIMA system that concentrations 
of chlorine dioxide, iodine and malonic acid vary little on the scale of variations in the 
chlorite and iodide concentrations ||. This observation has formed the basis of the adiabatic 
elimination of the former reactants, resulting in a two-variable reduction of the full model 



23] , and making their identification as the background species a natural one. Hence, in the 



interest of aligning experiment with numerical and theoretical work in this area, it appears 
reasonable for experiments to implement the CDIMA system. 



III. NUMERICAL METHOD 

A pseudospectral method was used to solve the governing partial differential equation in 
two dimensions. The physical boundary conditions are no-flux in the x-direction (transverse 
to the gradients), and fixed point boundary conditions in the z-direction (along the gradi- 
ents). In our numerical implementation, the governing equations are cosine Fourier trans- 
formed in the x-direction, and each spectral mode is evolved in time as a one- dimensional 
problem in the z-direction. We used a five-point finite-difference approximation with a 
variable-width spatial mesh to allow better resolution of the sharp front patterns along the 
gradients. The time-stepping scheme was Crank-Nicolson for the linear terms (implicit) and 
Adams-Bashford for the nonlinear terms (explicit), both second order accurate. After each 
time step, the updated solutions were transformed back into real space to reconstruct the 
nonlinear terms. 

As Table [I] indicates, the large order-of-magnitude variations in the real chemical pa- 
rameters of the LRE model make numerical solution of this model less tractable than the 
aforementioned abstract models. The stability restriction on the time step At due to explict 
treatment of the nonlinear terms is onerous. Hence, we parallelized the above numerical 
scheme, and the computations were performed on a 512-node Intel Paragon. 

As initial conditions, we used the stationary solution ut,(z) along z with a uniform profile 
in x, to which we added a perturbation given by the most linearly unstable mode, k*: 

~Tt(x, z, t = 0) = ~Tt s (z) + C8uk*{z) cos (k*x). (12) 

5uk*(z) is the eigenvector obtained from the linear stability analysis, and C is an overall 
scale factor to ensure that the perturbation is small and lies in the linear regime. The 
concentration vectors correspond to the six variables of the LRE model. The full six-variable 



linear stability analysis was carried out using inverse iteration |?1 . The non-zero boundary 
conditions in the ^-direction are: [MA] L = 0.0115 M, [C10 2 ] L = 0.006 M, and [I 2 ] L = 0.008 
M, where (R, L) refer the right and left reservoirs, respectively. From the linear stability 
analysis of ul(z), the growth rate for the instability at k* = 471.2 cm -1 is X(k*) = 0.00465 
s^ 1 , giving a characteristic saturation time of r ~ 1/A ~ 215 s. The system size is 0.3 cm in 
the z-direction and 0.133 cm in the x-direction, corresponding to exactly ten wavelengths 
in the x-direction. The spatial resolution of the system investigated here was N x = 129 and 
N z = 914. One-thousand iterations at this resolution required approximately 12 node-hours. 

The integration time step was At = 0.001 s, chosen to be the same as that for the time 
evolution in one dimension. In the one-dimensional time evolution, the restriction on the 
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time step due to the explicit treatment of the nonlinear terms was explored empirically, by 
varying At and using a value such that the algorithm was stable. We investigated using 
a higher order (third order) Adams-Bashford scheme to verify that the restriction on At 
was limited by stability as opposed to accuracy. The dynamically evolved stationary state 
was compared with that obtained from direct solution using a Newton-Raphson root-finding 
algorithm. We determined the time step used in one dimension to be adequate for the 
time evolution in two dimensions as follows: By using initial conditions uniform in the IE- 
direction (and equal to 5 x 10~ 12 M for all species), thereby reducing the two-dimensional 
time evolution to be effectively one-dimensional, we verified that the asymptotic solution was 
the same as that obtained in one dimension. It is possible that implementing a numerical 
algorithm adapted to solving stiff partial differential equations (see Ref. |12| and references 
therein) will improve the total execution time in the two-dimensional evolution. 



IV. RESULTS 
A. Two-dimensional stationary solution 

In Fig. [I], we show density plots of the initial state, as well as the numerical solution 
for the starch triiodide species after a total integration time of approximately T = 1000 s, 
with dark and light shading representing low and high concentrations, respectively. In Figs. 
|2](a)-(d), we show the time evolution of the k* mode and its higher order harmonics at the 
peak of the linear instability eigenvector (z = 0.094 cm) for the starch triiodide species. 
In Figs. §(e)-(h), we plot the logarithm of the magnitude of these quantities. (These plots 
are shown for the non-dimensionalized quantities.) The slope of the linear segment in Fig. 
|(e), for t < r ~ 0.2 is 5.129 ± 0.012 (corresponding to 0.004616 ± 0.000011 s) agrees well 
(to within one percent) with the linear growth rate, and further verifies our linear stability 
results. In these plots, it is apparent that the saturated amplitudes of the higher order 
harmonics are much smaller than that of k*. 

To compare the size of the nonlinear perturbations in the x-direction with the unper- 
turbed profile in the z-direction, we have plotted the x- and z-dependence of the two- 
dimensional solution. In Fig. 0(a), we show the x-dependence at the peak of the linear 
instability eigenvector (z = 0.094 cm), which approximates a pure Fourier mode, verified by 
the relatively small saturated amplitudes of the higher order harmonics in Fig. Q In Fig. 
H(b), the dashed and dotted lines denote the profiles in the z-direction at x — 0.067cm = 5A 
and x = 0.073cm = 5.5A, respectively. The solid line denotes the profile of the unperturbed 
one-dimensional stationary state. We note that although the saturated amplitude of the 
transverse instability is comparable to the variation of the one-dimensional stationary state 
in the z-direction, its x-dependence is not strongly nonlinear. 

The results presented in this section can be summarized by three points. First, we 
have presented the first numerical solution of two-dimensional patterns in the boundary- 
fed CDIMA reaction-diffusion system using the LRE model. Second, our numerical solution 
agrees qualitatively with patterns observed in thin-strip reactors in experiments on the CIMA 



system, for experimental conditions giving a single unstable front [25]. The wavelength of the 
solution presented here is 0.133 mm, in agreement with the 0.13 — 0.33 mm experimentally 
observed range f26j. Finally, the numerical evolution in two dimensions confirms our result 
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from the linear stability analysis of the one- dimensional stationary state along the gradients 
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B. Bifurcation diagram 

The symmetry-breaking transition from Fig. |l](a) to Fig. |I](b) is effectively one- 
dimensional, since only a single layer is unstable over a range of values of [MA]l control 
parameter, as we showed in reference In this earlier work, we demonstrated the exis- 



tence of (three) disjoint unstable ranges as the [MA]l control parameter was continuously 
varied from 4.0 x 10~ 3 M to 4.0 x 10~ 2 M, consistent with experimental observations in 
thin-strip reactors showing the appearance and subsequent vanishing of a transverse insta- 
bility as one of reservoir concentrations was increased. Here, we numerically investigate the 
dependence of the saturated amplitude of the transverse instability on the [MA]l control 
parameter in the vicinity of the (lower) bifurcation point for one of these unstable ranges 
(9.73 x 10~ 3 M < [MA]l < 1.25 x 10~ 2 M). In the following, for convenience, we report our 
results in non-dimensionalized units: the time conversion factor is k\ a = 9 x 10~ 4 s _1 , and 
the concentration conversion factor is k\b = 5 x 10~ 5 M. 

First, the critical control parameter was determined numerically from a linear fit to the 
maximum linear growth rate versus [MA]l- This value was found to be [MA] C = 194.5226. 
Linear stability analysis of the stationary state at this value of [MA]l yields a maximum 
growth rate of A* = 1.1054 x 10~ 4 . This value of A*, which is expected to be zero, gives a 
combined measure of the numerical uncertainties in the determination of the stationary state 
at a particular value of [MA]l as well its linear stability. Hence, in principle, there would 
be error bars associated with values of [MA]l, equal to A[MA]l = a AX = 3.0826 x 10~ 4 , 
where a is the slope of the linear fit from which [MA] C is extracted. 

Fig. |] shows the computed bifurcation diagram. Starting in the supercritical regime, 
for each value of [MA]l, as initial condition, we use the corresponding one-dimensional 
stationary solution in the ^-direction seeded with the most linearly unstable eigenvector at 
approximately the saturated amplitude of the previous higher value of control parameter 
(adiabatic approach). The final converged amplitudes were obtained from a least squares fit 
of the dynamical evolution of the peak amplitude to an exponential plus a constant offset, 
excluding initial points until the converged value did not depend on the number of excluded 
points. 

Close to onset, the solution is given by: 

u(x, z, t) = 5u k * (z) ■ [A(t) exp (ik*x) + A*(t) exp (-ik*x)] , (13) 

where k* is the most linearly unstable mode, and A(t), A*(t) are complex conjugates. The 
amplitude equation, which gives a universal description of the weakly nonlinear behavior 
and depends only on the symmetries of the problem (in this case, translational invariance 
in the ^-direction) is (at seventh order): 

8 A 

— = eA + 9l \A\ 2 A + g 2 \ A\ A A + g 3 \A\ 6 A, (14) 

where the coefficients depend on the specific system under investigation, and e = [MA]l — 
[MA] C is the distance from onset. Our results, described below, reveal a subcritical (first 
order) bifurcation. 
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Fig. £| shows a sixth order polynomial fit to the numerical data: [MA]l = [MA] C — g\ \A\f — 
g 2 \A\ 4 -g 3 \A\ b . First, [MA] C was held fixed at the linear threshold value, and (gx, g 2 , g 3 ) were 
fitted for. The goodness of the fit depended on the number of points farthest from linear 
threshold retained in the fit. We determined that excluding the last point ([MA]l = 230.0) 
gave the best fit. The fitted parameters are: 

gi = 1.7585, (15) 
g 2 = -4.0825, (16) 
g 3 = -0.82677. (17) 

Since g\ > 0, the instability does not saturate for e > 0, and the bifurcation is subcrit- 
ical. We also investigated allowing all parameters, ([MA] C , gi, g 2 , #3), to float. Again, ex- 
cluding the last point produces a fit with an offset [MA] C , which is closest to the linear 
threshold value. The values of the fitted parameters in this case are: ([MA] C , g\, g 2 , #3) = 
(194.47, 1.5548, -3.8692, -0.89160). 

Fig. [5] shows the convergence of the peak amplitude of the fastest growing linear mode 
for [MA]l = 194.4 in the subcritical region. It shows convergence from above and below to a 
finite amplitude instability. The error for the converged amplitude is estimated to be one-half 
of the difference between the converged-from-above and converged-from-below values. This 
error (8.0 x 10 -4 ) is taken to be the same for all points, eventhough the convergence from 
below was not repeated for all pointsfj We also confirmed the decay of a linear perturbation 
at this same subcritical value. An exponential fit to the dynamical evolution of the peak 
amplitude of the instability yields a decay rate of A = —0.0416, in good agreement with the 
largest eigenvalue, A = —0.0429. 

The minimum value of [MA]l below which a finite amplitude instability does not exist, 
corresponding to the saddle-node bifurcation, can be computed from the fitted parameters. 
Using the parameter values given in Eq. p~5| [TT|, [MA]sn is found to be: 

[MA] SN = 194.34. (18) 

The inset in Fig. ^ shows this turning point. For [MA]l = 193.0 below this value, we 
explicitly verified decay to zero of an initial perturbation with amplitude A = 0.7756. 

We note that the transition is "weakly" subcritical. This is characterized by the small 
range of control parameter below linear threshold, approximately equal to 0.18 (9 x 10~ 6 M), 
for which a finite amplitude instability exists, in comparison with the linearly unstable range, 
55.4 (2.77 x 10" 3 M), determined in our earlier work [O]. The weakly subcritical nature of 



the transition implies that a linear stability analysis of the one-dimensional structures along 



the gradients flTIj] does have utility in predicting the existence of a transverse instability for 
given reaction parameters and boundary conditions over a wide range of control parameters 
(in the supercritical regime). 



1 This convergence from above and below was also confirmed for a point in the supercritical region, 
[MA]l = 195.0, with even better agreement between the two fitted converg ed values (1.0 x 10~ 4 ). 
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V. CONCLUSION AND DISCUSSION 



To summarize, we have carried out a two-dimensional numerical simulation of the 
boundary-fed CDIMA react ion- diffusion system based on the realistic LRE model for this 
system. Our results are qualitatively similar to those seen in experiments on the CIMA 
react ion- diffusion system, and support the earlier work |TT[ in which we studied the linear 
instability of the boundary-fed CDIMA system to transverse Turing patterns. 

Numerical studies by Jensen et al. p7|-p9[ on the two-variable LRE model with uniform 
backgrounds have found the transition to stripes in one and two dimensions to be sub- 
critical. Our results demonstrate this transition to also be subcritical in the boundary-fed 
LRE model. This prediction and computed bifurcation diagram can be directly verified by 
experiments based on the CDIMA system. 

The subcritical nature of the transition to stripes makes the LRE model qualitatively 
different from other abstract reaction-diffusion models hitherto used to study Turing pat- 
terns. For example, in the ramped Brusselator, the transition to stripes has been shown to 
be supercritical |3(| . Jensen et al. have investigated the propagation of fronts separating the 
homogeneous steady state from the Turing structure in one and two dimensions using the 
uniform LRE model. The subcriticality allows for the existence of a range values of control 
parameter for which the front velocity vanishes, allowing an infinite number of stable steady 
inhomogeneous structures. Despite the weakly subcritical nature of the transition, it would 
be interesting to similarly investigate front propagation and formation of localized (quasi 
one-dimensional) states in the boundary-fed system. 

In experimental geometries (disc reactors) where the dimensions of the reactor transverse 
to the gradients are large, the analog of the one- dimensional row of spots which develops in 
our numerical simulation and in experiments using thin-strip reactors is a two-dimensional 
"monolayer". Dufiet et al. |JT] have pointed out that these monolayers, which are confined 
by a strong transverse gradient of reservoir chemical concentrations, must be distinguished 
from genuine two-dimensional structures with uniform control parameters. Pattern selection 
in genuine two- and three-dimensional systems has been studied analytically and numeri- 



cally using abstract reaction-diffusion models P%fl . However, it is not practical to generate 
sustained genuine structures experimentally. In the context of a model reaction-diffusion 
system, Dufiet et al. have shown that in genuine two-dimensional systems and monolay- 
ers, the stripe-hexagon competition is similar close to onset. They find, however, that far 
from onset, hexagonal phases in monolayers are restabilized due to their interaction with 
a longitudinal (k = 0) instability. The latter finding is consistent with earlier theoretical 
predictions p3| , |3^| , as well as experiments in bevelled disc reactors |35 |. 

It would be interesting to numerically investigate pattern selection for monolayers in the 
LRE model of the CDIMA system in the range of boundary conditions and reaction param- 
eters accessible to experiments, allowing in principle direct comparison with experimental 
results. This would require extension of our numerical computation to three dimensions. 
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TABLE I. Kinetic constants for the CDIMA system. 



Rate or diffusion constant 


Dimensions 


Value 


ha 


(s- 1 ) 


9 x 10" 4 1 


hb 


(M) 


5 x 1(T 5 1 


h 


(M -1 s -1 ) 


1 x 10 3 1 


ha 


(M-V 1 ) 


1.2 x 10 2 1 


hb 


^ (s- 1 ) 


1.5 x 1(T 4 1 


h 


(M 2 ) 


1.0 x 10~ 14 1 


k + 


(M- 2 s~ 1 ) 


6.0 x 10 5 2 


k- 


(s- 1 ) 


1.0 2 


D\- 


(cm 2 s _1 ) 


7.0 x 10~ 6 3 


D cio 2 - 


(cm s J 


7.0 x 10~ 6 3 


D h 


(cm 2 s _1 ) 


6.0 x 10~ 6 1 


Dma 


(cm 2 s _1 ) 


4.0 x 10~ 6 1 


Dc\o 2 


(cm 2 s _1 ) 


7.5 x 10~ 6 1 


Dm 


(cm 2 s _1 ) 


1.0 x 10~ 5 


K[S} 


(M- 1 ) 


6.25 x 10 4 4 


Trom (I| at 7°C; 2 From f2 


at 4°C; 3 From [|| at 4°C; 4 From 


at 4°C. 
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z (cm) 





0.075 



0.15 



0.225 



0.3 



z (cm) 





0.075 



. 15 



0.225 



0.3 



FIG. 1. Two-dimensional solution for starch triiodide with [MA]l = 0.0115 M: The top figure 
shows the numerical solution after T = 1000 s of evolution time. The bottom figure shows the 
initial condition. Dark and light shadings correspond to low and high concentrations, respectively. 



13 




0.5 1.0 1.5 



FIG. 2. Time-evolution of the most linearly unstable mode k* and its higher order harmonics: 
(a)-(d) show the values of the k*, 2k*, 3k* and 4k* modes for the starch triiodide species at the 
peak of the most linearly unstable eigenvector (z = 0.094 cm), as a function of time; (e)-(h) show 
the logarithms of magnitudes of these amplitudes. The plots are shown for non-dimensionalized 
quantities. (The time and concentration conversion factors are 9 x 10~ 4 s" 1 and 5 x 10~ 5 M, 
respectively.) The slope of the linear segment (heavy line) in plot (e) for t < 0.2 is 5.129 ± 0.012, 
and agrees well with the growth rate X(k*) = 5.172 from the linear stability analysis. 
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FIG. 3. Spatial dependence of the two-dimensional solution in the x- and z-directions for 
SI,3~ species: (a) The system size in the x-direction is 0.133 mm, constructed to incorporate 
exactly ten wavelengths of the most linearly unstable mode k*. The profile in the x-direction is 
plotted at z = 0.094 cm. (b) The dashed and dotted lines denote the profiles in the z-direction 
at x = 0.067cm = 5A and x = 0.073cm = 5.5A, respectively. The solid line denotes the profile of 
the unperturbed one-dimensional stationary state. We note that the saturated amplitude of the 
instability is comparable to the variation in the profile of the one-dimensional stationary state in 
the z-direction. 
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FIG. 4. Bifurcation diagram: The solid curve is the computed fit; the broken curve cor- 
responds to the unstable branch. The inset shows the vicinity of the saddle node bifurcation 
([MA] S n = 194.34), and the dotted vertical line denotes linear threshold ([MA] C = 194.5229). 
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FIG. 5. Convergence to finite amplitude below linear threshold, [MA]l = 194.4: The 
closely-spaced circles denote the numerical time evolution, and the solid lines denote the com- 
puted fit to an exponential plus a constant offset. Convergence from above and below to a finite 
amplitude is apparent. 
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